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Abstract 

The objective of this paper is to prove the convergence of a linear implicit multi- 
step numerical method for ordinary differential equations. The algorithm is obtained via 
Taylor approximations. The convergence is proved following the Dahlquist theory. As an 
additional topic, the time stability is established too. Comparative tests between some of 
the most known numerical methods and this method are presented. 

Keywords Taylor approximation, multi-step method, stability, consistency, time regions sta- 
bility. 
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1 Introduction 

We present a linear implicit m-step method LIL (Local Iterative Linearization) and prove 
its convergence applied for the following initial value problem 

x^f{t,x), x{to)^xo, (1) 

where / : [io, T] x M" ^ E", T > 0, io S K+, is a C" smooth Lipschitz functiorQ 

Although the classical linear multi-step algorithms are very known and utilized, the LIL 
characteristics (convergence properties, time stability and applications results) show that this 
numerical method could be considered as an interesting alternative to the widely used formulas. 

The backward approximation of derivatives implies null coefficients of the odd order deriva- 
tives which represent a major advantage for the propagation of errors. 

As a comparative test two simple ODEs with known analytical solutions and a chaotic 
continuous-time dynamical system, first studied by Fabrikant and Rabinovich [6] and recent 
numerically re-examined by Danca and Chen [3], was integrated using the LIL algorithm and 
some of the most known algorithms. The complex dynamic of this special model represented 
a real challenge for almost all of these methods as shown in Sect. 5. 

Being an implicit method, an extrapolation is used as the predictor phase. Like all the m- 
step algorithms, the previous m points (beside the first m start points) should be estimated 
every step. 

^The Lipschitz condition is necessary for the stability proof. 
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To study the convergence we use the unified approach of stability and consistency developed 
by Germund Dahlquist in 1956 [2] (see also [7-8]). Thus, the LIL method applied to the initial 
value problem ([T]) is considered convergent if and only if it is stable and consistent. 

The content of this paper is as follows: In Sect. 2 the LIL method is deduced. The conver- 
gence is proved in Sect. 3. In Sect. 4 is presented the time stability with the corresponding 
time stability domains. In Sect. 5 three examples are presented. All computer tests were 
realized using a Turbo Pascal code written by the author. In the Appendix, the coefficients of 
LIL method are presented. 

2 Deduction of the LIL method 



Let us consider the uniform grid 

A = (to < ti < ... < t„ = T), neN* 

with the step-size 



T-tn 



2 St. 



where 6t stands for the ray of the neighborhood Vk = {tk — St , tk + St) , k — 1, 2, n — 1. 
We assume that all infinite Taylor series converge, but this is not necessarily since one truncate 
at a sufficiently large but finite number of terms. 
One introduce the following notations 



= x(tk-j) ^x[to + {k~ j)h)] , 

= x«(t,), xi°^ ■.= x{tk), J = 1,2,... 



In the following k is supposed to take the values k — 1,2, ...,n — 1. 

If we consider function of variable h defined in Vk , then the first m terms of 

Taylor approximation of Xk^j i^ 



Xk- 



Xk 



2! 



(-irl^4™\ 



where j = l,2, ...m. The relations (^2| represent a Cramer system with the unknown x 



(2) 



i) ■ 

fe ' « = 



1, 2, m 



Xk~l Xk ~ i\^k~^ 2! -^fe ••• + ( 1) m{^k ^ 



Xk-2 ~ ^k 



2h' 
1! -^k 



(2/t)% " / (2/t)"' (m) 

2! -^k ■■■ \ ^) m\ •''k ' 



(3) 



Xk—m Xk 

The determinant of the system ([S]) is 



nA' I (™^') J' 
1! -^fc 2! k 



A 



■^(n^+l)/2 



1!2!...to! 



1 1 1 

2 22 23 



1 

2" 



+ (-1) 



m (m h)'^ (m) 
m ! k 



m(m+l)/2 

= h l!2!...(m-2)!, 



^The choice of m and h is supposed to be such that the Taylor approximation can be used. The link 
between h and m is analyzed in Section 3.1 
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Because for m > 2 we have A 7^ 0, there exists a unique solution 

1 
h 

for m = i = 1 



^'^ - ^"^SiyXk-j, i = 2,...,m for m>i>l, 



Xk " Xk-l 

h 



the coefficients Sij being drawn in Table 7/ Appendix. 

Thus we obtained a backward approximation of derivatives, which represents the key of LIL 
method. The Taylor approximation of the solution x, considered now as function of t in the 
neighborhood Vk , is 



1! 2! m! 



(4) 



Next, integrating Q in Vk we get 

tk+St St St 

I ■ 

tk — S t 



k+St St / ™ r> \ 

/ x{t)dt^ J x{t + tk)dt» J E 4 

.-St -St -St\i=0 J 



1=0,2,4, 
i< m 



V 1 ui+lji)^u„, I /)3 1 „" I u5 1 ^(4) I 

^ 2'{i+iy."' ■''k ^llXk-tlL 24'''fe 1920'''fe ^ 



(5) 



Remark 1 T/ie zero coefficients of the derivatives x^k^'^^'' for i — 0,1,2,... in ^ represent 
a major advantage for the propagation of errors and computation time. 

If we use in ^ the derivatives expression (2.3) we have 

tk+St 
^ p m 

/ x{t)dt » h'^aoiXk-i (6) 

tk — o t 

the coefficients a^i being given in Table 8(a)/Appendix. Using the same way one can ap- 
proximate x' on Vk 

tk+St 

/ x{t)dt»'^aiiXk-i, (7) 
tk'St »=o 
the coefficients an being drawn in Table 8(b) /Appendix. 

To overcome the difficulty of Taylor approximation of the composite function / we found, 
empirically, that the relations ^ could be considered as a simple way to approximate the 
integral of / without altering the method convergence. Thus 

J f{t,x{t))dt^hY,<^O^fk-^, (8) 

tk-St ' = 

where fk-i ■= f {tk-i,x{tk-i)) ■ 

Using ([7]) and (|8| we can integrate ([T]) in Vk 

m m 

cTiiXk-i = h ^ aoi fk-i ■ 

i=0 i=0 
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m 




1 


Xk = Xk-i + hfk, 


2 


Xk = ^Xk-i - |xfe„2 + 5|(25 /fc - 2 fk-i + fk-2) , 


3 


Xk = ^Xk-i ~ ifa;fc-2 + lxk-3 + j| (26 A- - 5 A-i +4 A-2 - A-a) 


4 


Xk = 2xk-i - |xfe_2 + §Xk-3 - |2;fc-4 + Y^{Q463 fk - 2092 A-i 
+2298A-2 - 1132A-3 + 223A_4), 


5 


Xk = Ixk-l llxk-2 + llxk-3 tl^k-A + lxk-5 + 14^75 (6669 A 

-3122 A-1 + 4358 A-2 - 3192 A-3 + 1253 A-4 - 206 A-s)- 



Table 1: LIL algorithms. 



Because aiQ ^ , for every m (see Table 8/ Appendix), the approximation of the solution 
in Vk is 

h 

Xk 



If we denote 



o-Q i fk-i — i xk-i ■ (9) 

2 — 2 — 1 

ECTOi A-ii ''^fc / , ' 



CTin CTin 

^^^2 = ^'^ 2=1 



the relations ([9| become 

Xk^Vk + huk, fc = 1, 2, ...n - 1 . (10) 



Formula (10 1 represents the mth-order LIL method. In Table 1 the formulae for orders one 
through five (m G {1, 2, 3, 4, 5}) are presented. 

The study was achieved up to m = 8, but in this paper for the sake of simplicity we 
considered only m e {1, 2, 3, 4, 5}. For m 1 the LIL method is equivalent to the backward 
Euler method. 

The LIL method is an implicit method due to the presence of the term A in the right 
hand side which depends on Xk- Therefore additional computations are necessary in order to 
calculate A- In this purpose we approximate Xk-i £ Vk (see ([2])) 



Xk-i xxk- Y^Xk + -^Xk - ... + (-1) —^xl '. 
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Using for derivatives the relations (2.3) one obtains 



wherefrom we have 



2 i-^k—i 



i=0 



i=0 



m i-^k—i 5 



Xk 



'^EmiXk-i, « = 1,2,...,TO, m > 1. 



(11) 



The coefficients Smi are given in Table 9/ Appendix. 
Using (111, fk becomes 

/ m ^ 



The relation (111 represents an extrapolation formula (predictor phase) for Xk and can 
be used to approximate the solution, but without an acceptable accuracy, while ( 10 1 is the 
corrector phase. 

Because ( 10 1 is a multi-step relation, a starting method (for example the standard Runge-Kutta 
method) is necessary in order to calculate the m first start values: a;_i,a;_2, ...jX-m- 



3 The convergence 

The convergence is analyzed using the Dahlquist theory which states that a numerical 
method is convergenlj^ if it is consistent and stable (see [2], [4] or [7-8]). In this purpose let us 
consider the LIL method ( 10 1 in the usual form 

<7wXk + CrilXk-l + ■■■ + CTimXk-rn = CToo fk + CTOl fk-1 + ■■■ + (^Om fk-m, (12) 

with the characteristic polynomials 

m m 

=E'^i'^'""'' /?™(.s) = 5]ao,s™-\ me{l,2,3,4,5}. (13) 

3.1 Consistency and errors 

Following the Dahlquist theory, the LIL method is consistent because its characteristic 
polynomials (131 satisfy Q!m(l) — and a^(l) — — /3„i (1) for m e {1,2,3,4,5}. As it 
is known, the order of a linear multi-step method is r if, and only if, r of the following 
coefficients 

m m 

vanish. 

Note that above the convention O" = 1 was used. The values of C for LIL method are given 
in Table 2. 

From Table 2 one can deduce that the LIL order (the largest r for which C is null) is m + 1. 
The local truncation error ej is, for a given to, of order to -I- 1 (see e.g. [7]). 



''The "convergence" means here "uniform convergence" on an interval for any C" smooth function /. 
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m 


Ci 


C2 


C3 


Ca 




^6 






1 








-0.5 












2 











-0.04 










3 














-0.313 






0{h^) 


4 

















-1.37 




Oih'') 


5 




















-177.184 


0{h'') 



Table 2: C coefficients. 



m = 2 


7i(s) = (3s-l), 


m = 3 


72(s) = (15 - 10s 3), 


m = 4 


73(s) = (35 - 35 21 s - 5), 


m = 5 


74(s) = (315 s-* - 420 s'^ + 378 - 180 s 35). 



Table 3: The polynomials jm-i- 



Comparatively, the local truncation error for the standard (4th-order) Runge-Kutta algo- 
rithm is of order 4, and for the multi-step algorithms Adams-Moulton and Gear are of order 
m + 1, the same as for LIL algorithm. 

The global truncation error (the accumulation of the local truncation errors) per unit time 
is e^" = £t/h. Hence the global truncation error per unit time is of m order. 

3.2 Stability 

LIL is stable if all solutions of the following difference equations 

a,„(s) = 0, me {1,2,3,4,5}, (14) 

are bounded. A necessary and sufficient condition for stability is that all zeros Sfe , fc = 
l,2,...,m of am satisfy |sfc| < 1 and that zeros with |sfc| = 1 be simple. It is easy to 
see that ai{s) = s — 1 and for m > 2, am{s) = (s — l)7,„_i(s) (Table 3) with the zeros, 
numerically found for m = 3, 4, 5, given in Table 4. 

Hence the LIL method is stable and therefore we have the following result 

Theorem 2 The LIL method for to the initial value problem ([T]) is convergent for all m G 
{1,2,3,4,5}.. 

Proof. Because LIL is consistent and stable, following the Dahlquist theory, it is conver- 
gent. 





Sl 


S2 


S3 


S4 


S5 


m = 2 


1 


0.33 








m = 3 


1 


0.33 -hi 0.30 


0.33- i 0.30 






m — A 


1 


0.40 


0.30 + i 0.52 


0.30 -i 0.52 




m = 5 


1 


0.40 + i 0.17 


0.40 0.17 


0.26 -hz 0.72 


0.26 -i 0.72 



Table 4: The zeros of the characteristic equation (3.3). 
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4 The regions of time stability 



An integration method may have low round-off error and low truncation error, but be totally 
worthless because it is time unstable. The standard method for testing the time (numerical) 
stability is to apply the integration method to the first-order linear test equation 

X = Ax, x{Q) = xo, (15) 

where x, a;o,A may be complex. A method is time (numerically) stable for specified values 
(A, h) if it produces a bounded sequence {a;„} when applied to the test problem ( [T5| [7]. The 
set of the complex values z = Xh for which {xn} is bounded is called the stability region of 



the method. When an integration method is applied to the system ( 15 ) the result is a linear, 
discrete-time system with a fixed point at the origin. This means that the stability regions 
contain the half plan Re{z) < 0. Therefore the stability of this fixed point determines the time 
stability of the integration method. 

Although this stability criterion guarantees that a method is stable only when integrating a lin- 
ear system, and not for nonlinear systems it is an usual way to compare numerical performances 
for different algorithms. 

Following the theorem which states that a linear multi-step method is time stable for a 
particular z if and only if, the equation a„i(^) = z has the following properties: all 
roots satisfy | f | < 1, and all roots with | ^ | = 1 are simple (see e.g. [8]), the proof of the 
time stability of LIL method follows from convergence study. 

In order to draw the stability regions let us define 

PrniO :-a„,(e)-z/3™(e), 

Then, a linear multi-step method has the stability region S , the set of all points z G C 
such that all the roots of Pm{£,) — lie inside or on the unit circle and those on the unit 
circle are simple. Hence we obtain the equation 



which has to be solved for any given z d C . But instead of solving (16) for given z , we can 
give ^ = e*^ with | ^| = 1 and plot 

am(e'^) , . 

/3,n(e»'')' ^ ' 

for e [0, 27r] The set thus mapped must contain dS . The stability region of a numerical 
stable algorithm has to contain the origin in his boundary. 

In Figure 1 the stability regions for LIL algorithm for m S {1,2,3,4,5} are drawn. One can 
observe that LIL algorithm has, for all m, large (even unlimited) regions of stability, including 
the entire left-half complex plane, typically for implicit algorithms. The time stability of LIL 
method is more efficient than that of other known algorithms and is comparable with time 
stability of the Gear's algorithm (see e.g. [5] where the stability regions were drawn for several 
known algorithms). 

Taking account of the fact that higher order is not always higher accuracy, an acceptable 
compromise between the accuracy, time stability and computational time was proved to be 
TO = 3. 
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5 Applications 



5.1 LIL versus standard methods 

The goal of this section is to compare the characteristics of few known standard algorithms (the 
4th-order methods: Runge-Kutta, Gear, Adams-Moulton, the 3th-order Adams-Bashforth 
method and the Milne method) and 4th-order LIL method. For this purpose we integrated 
two simple examples, with known analytical solutions: the Bernoulli equation 

2fx{t) - Atxit) - x^{t) = 0, 

and 

X [t) — cos(i). 

The following values were calculated: 

- the relative error £r = — x\ / "Y^Xa , where Xa is the analytical solution. The sum 
is taken over the integration interval. 

- the maximum absolute error: A = max|a;a,A; — Xk\, where Xa,k is the exact solution in 

- the computation time ^ 

The results are presented in Table 5 and 6. 





(a) 


R-K 


Gear 


A-M 


A-B 


Milne 


LIL 






L9-10-'' 


LQ-IO-'^ 


1.1-10-'' 


1.1-10-^ 


1.4-10-"^ 


1.4-10-^ 


A 


L9-10"2 


2.0-10-2 


1.2-10-^ 


1.2-10-^ 


1.8-10-^ 


1.5-10-'^ 




0.16 


0.16 


0.16 


0.10 


0.10 


0.16 


(b) 


R-K 


Gear 


A-M 


A-B 


Milne 


LIL 




3.8-10-^ 


3.8-10-^ 


2.3-10-1" 


2.3-10-1" 


2.8-10-1" 


2.8-10-1" 


A 


1.9-10"^ 


2.0-10"^ 


1.2-10-** 


1.8-10-« 


1.8-10-** 


1.5-10-** 




0.82 


0.71 


0.82 


0.43 


0.71 


0.87 



Table 5: Bernoulh equation integrated with: a) /i = 0.01, t e [1, 100]; b) /i = 0.001, t e [1, 50]. 



Comparing the results in Tables 5 and 6 one can deduce that LIL's performances, for these 
two examples, are comparable to those of performant methods like Gear, Adams-Moulton and 
Adams-Bashforth. 

5.2 Rabinovich-Fabrikant system 

The hard test was the integration of the Rabinovich-Fabrikant system. Rabinovich and 
Fabrikant [6] studied the following dynamical system (named the R-F model hereafter) 

2^1 — X2{X3 — 1 + xf) + axi, 

X2 = xi{3x3 + 1 - xl) + ax2, a, 6eIR. (18) 

X3 = -2X3(6 + X1X2), 

is here only a relative value since it depends on the used code (Turbo Pascal and using 64 bits), and the 
computer processor (500 MHz ). 



8 



(a) 


R-K 


Gear 


A-M 


A-B 


Milne 


LIL 


£r 


2.4-10-^ 


4.9-10-^ 


7.7-10-'^ 


4.0-10-^ 


5.0-10-^ 


5.0-10-^ 


A 


3.7-10-^ 


5.3-10-^ 


4.9-10-'' 


2.6-10-^ 


3.9-10-^ 


3.3-10-^ 


t[s] 


0.0 


0.0 


0.0 


0.0 


0.0 


0.0 





R-K 


Gear 


A-M 


A-B 


Milne 


LIL 




4.9 - lO-'' 


9.9 - 10-4 


3.9 - 10-7 


1.5-10-" 


1.9-10-** 


1.9 - 10-** 


A 


7.5 - 10-4 


1.0 - 10-^ 


2.4-10-7 


2.7-10-7 


1.5-10-** 


1.2 - 10-6 


t[fi] 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 



Table 6: x{t) = cos(i), t G [0, 27r] integrated with: a) /i = 0.05; h) h ^ 0.001. 



This system models the stochasticity arising from the modulation instability in a non-equilibrium 
dissipative medium. Some qualitative analysis and numerical dynamics have been reported in 
[6] and a carefully re-examination together with many new and rich complex dynamics of the 
model, that were mostly not reported before, are presented in [3]. The chaotic R-F model 
proved to be a great challenge to the classical numerical methods, most of them being not 
successful to study the complex dynamics of this special model. 

All computer test results and graphical plots in Figures 2-5 were obtained with a special Turbo 
Pascal code which plots phase diagrams and time series The code for LIL method may be 
obtained directly from the author. 

For a < b, the system is characterized by the appearance of chaotic attractors in the phase 
space (see e.g. Figures 2). 

It is well known that because of the sensitive dependence on initial data, a chaotic system tends 
to amplify, often exponentially, tiny initial errors. These kind of errors could be amplified to 
so large, that it is almost impossible to draw mathematically rigorous conclusions based on 
numerical simulations. A typical case can be seen from Figure 3, wherefrom one deduces that 
the attractor's size along the Xa-axis increases significantly as the step-size decreases. This 
problem has been noticed for a long time, and has promoted a useful theory called "shadowing," 
namely, the existence of a true orbit nearby a numerically computed approximate orbit [1]. 
We have also found that the strong dependence on the step-size for R-F system, for certain 
values of b and with the same initial conditions, could produce totally different attractors (see 
Figures 4) 

There are few special cases which proved to be a real challenge for the numerical methods. 
As example for the case a = 0.3 and b = 0.1 (shown in Figure 3), the 4th-order Runge- 
Kutta and Milne methods failed while only the Gear and Adams-Moulton methods seem to 
give comparable results to those obtained with LIL method; the attractors obtained with the 
3th-order Adams-Bashforth method are different to those obtained with Gear, Adams-Moulton 
and LIL methods (Figures 5). 



6 Concluding remarks 

In this paper we present a linear implicit multi-step method, LIL, for ODEs proving its 
convergence, too. The method could be considered as an acceptable alternative to the classical 
algorithms for ODEs and can be successfully used in practical applications. One of the advan- 
tages is that in ([5| only the even order derivatives appear, this fact reducing the truncation 
error and the computational time. 
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The algorithm seems to be stiffly-stable since it can integrate efficiently and accurately 
enough dynamical systems like R-F which presents stiflF characteristics. 

The implementation of adaptive step-size represents a task for a future work. The basic 
approach would be applicable directly to variable step-size. 
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(a) 




Figure 2: Two chaotic trajectories of R-F system: a) Three-dimensional phase portrait for 
a = 0.1, b = 0.2876; b) Plane phase portraits and time series for a = —1, b = —0.1. 
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(a) 




Figure 4: Two different attractors (plotted here by points), with the same initial conditions 
and parameters values (a = 0.12, b = 0.05), but with different step-size a) h = 0.05 and b) 
h = 0.005. 
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Figure 5: The case a = 0.3 and b = 0.1 integrated with: a) the Ath-order Adams-Moulton ; 
b) the 4i/i-order Gear algorithm; c) 3th Adams-Bashforth algorithm; d) 4th-LlL algorithm. 
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Appendix 



i = 2 


S20 


(521 


^22 


(523 


52A 


52b 


m = 2 


1 


-2 


1 








m = 3 


2 


-5 


4 


-1 






m = 4 


35/12 


-26/3 


19/2 


-14/3 


11/12 




m = 5 


15/4 


-77/6 


107/6 


-13 


61/12 


-5/6 



i = 3 


^30 


^31 


h2 


(^33 


(534 


<535 


m = 3 


1 


-3 


3 


-1 






m = 4 


5/2 


-9 


12 


-7 


3/2 




m = 5 


17/4 


-71/4 


59/2 


-49/2 


41/4 


-7/4 





i = 4 


(540 


5ai 


<542 


1^43 


(544 


(545 




m = 4 


1 


-4 


6 


-4 


1 




= 3 




-ii 


2() 


-21 


ii 


-2 




i = 5 


(550 


(551 


<552 


^53 


<554 


(555 




m = 5 


1 


-5 


10 


-10 


5 


-1 



Table 7: ^ coefficients. 



(«) 


TO = 1 


TO =2 


TO =3 


m =4 


m = 5 


foo 


1 


25/24 


13/12 


6463/5760 


741/640 


CTOl 





-1/12 


-5/24 


-523/1440 


-1561/2880 


(^02 




1/24 


1/6 


383/960 


2179/2880 


'^'03 






-1/24 


-283/1440 


-133/240 


^04 








223/5760 


1253/5760 












-103/2880 



(6) 


= i 


;/( =2 


;// =:! 


;/) =1 


=5 


CTlo 


1 


3/2 


15/8 


35/16 


315/128 


CTii 


-1 


-2 


-25/8 


-35/8 


-735/128 


0'12 




1/2 


13/8 


7/2 


399/64 


0-13 






-3/8 


-13/8 


-279/64 


ai4 








5/16 


215/128 












-35/128 



Table 8: a coefficients. 



TO 




£m2 


^7nZ 


£m4 




2 


2 


-1 








3 


3 


-3 


1 






4 


4 


-6 


4 


-1 




5 


5 


-10 


10 


-5 


1 



Table 9: e coefficients. 
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